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The gravitational radiation emitted dnring the merger of a black hole with a neutron star is 
rather similar to the radiation from the merger of two black holes when the nentron star is not 
tidally disrupted. When tidal disruption occurs, gravitational waveforms can be broadly classified 
in two groups, depending on the spatial extent of the disrupted material. Extending previous work 
by some of us, here we present a phenomenological model for the gravitational waveform amplitude in 
the frequency domain encompassing the three possible outcomes of the merger: no tidal disruption, 
“mild” and “strong” tidal disruption. The model is calibrated to 134 general-relativistic numerical 
simulations of binaries where the black hole spin is either aligned or antialigned with the orbital 
angular momentum. All simulations were produced using the SACRA code and piecewise polytropic 
neutron star equations of state. The present model can be used to determine when black-hole binary 
waveforms are sufficient for gravitational-wave detection, to extract information on the equation of 
state from future gravitational-wave observations, to obtain more accurate estimates of black hole- 
neutron star merger event rates, and to determine the conditions under which these systems are 
plausible candidates as central engines of gamma-ray bursts, macronovae and kilonovae. 

PACS numbers: 04.25.dk, 97.60.Jd, 97.60.Lf, 04.30.-w 


I. INTRODUCTION 

The year 2015 will mark the beginning of the advanced 
gravitational-wave (GW) detector era. Exactly one hun¬ 
dred years after Einstein formulated the theory of Gen¬ 
eral Relativity, the two Advanced Laser Interferometer 
Gravitational-Wave Observatory (LIGO) detectors [1, 2] 
are about to start their observation runs. They will soon 
be followed by Virgo [3], and later on by the Kamioka 
Gravitational wave detector (KAGRA) [4, 5] and LIGO- 
India [6]. Detections will provide us with unprecedented 
information about astrophysical GW sources. Coalescing 
compact binary systems containing neutron stars (NSs) 
and/or black holes (BHs) are the main target for GW 
interferometric detectors. Their waveforms encode in¬ 
formation about the masses, spins, distance, sky loca¬ 
tion, and orientation of the source, and, when NSs are 
present, about the NS equation of state (EOS). Detect¬ 
ing GWs emitted by compact binaries relies on match¬ 
ing noisy detector data with theoretical signal predic¬ 
tions and, therefore, requires us to build waveforms for 
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the targeted sources that are as accurate as possible. 
At the same time, interpreting future observations calls 
for understanding as many details as possible about the 
sources. 

In light of this, numerical relativity has made giant 
steps forward over the last decade, and simulations of 
the late inspiral and merger of compact binaries are now 
possible. As these calculations are resource intensive and 
time consuming, simulations that cover as many cycles 
as are necessary to fill the sensitivity bandwidth of the 
detectors and that span the whole parameter space are 
still beyond the reach of present-day computers. This 
is why semianalytical waveform models that hll the gap 
between perturbative methods — that describe the early 
inspiral stage — and numerical relativity are necessary. 
These models and numerical simulations are most ad¬ 
vanced for BH-BH systems. Phenomenological inspiral- 
merger-ringdown (IMR) waveform models have been pro¬ 
posed for nonspinning binaries by Ajith et al. [7-9]; for 
spinning, nonprecessing binaries by Ajith et al. [10] and 
by Santaman'a, Ohme, et al. [11], and more recently by 
Khan et al. [12, 13]; and for spinning, processing binaries 
by Hannam et al. [14]. These are generally referred to 
as “PhenomA,” “PhenomB,” “PhenomC,” “PhenomD,” 
and “PhenomP,” respectively, and are all based on a 
post-Newtonian (PN) description of the early inspiral. 
Similarly, a lot of effort was put into tuning the effective- 
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one-body (EOB) model to BH-BH numerical simulations, 
and progress in this direction has been remarkable (see 
e.g. [15-21]). IMR models are employed in several con¬ 
texts: to design and build GW detection templates at 
manageable computational costs (e.g. [22, 23]), to test 
GW detection infrastructure (e.g. [24-26]), to evaluate 
statistical and systematic errors (e.g. [27, 28]) and to per¬ 
form detection rate calculations (e.g. [29-32]). 

In the case of compact binaries containing at least one 
NS, the waveform modelling landscape is less developed, 
both because long and accurate simulations are particu¬ 
larly hard to achieve and because the parameter space is 
larger. The outcome and aftermath of NS-NS and NS-BH 
binary simulations (as opposed to BH-BH binaries) de¬ 
pends on several assumptions on physics that is currently 
underconstrained, including e.g. the NS EOS, the effect 
of magnetic fields and neutrino emission. This makes the 
NS-NS case particularly complex, as a hypermassive NS 
may form in the merger and oscillate for 10-100 millisec¬ 
onds, emitting GWs in a fashion that is hard to predict 
from the parameters of the binary itself [33-35]. Due 
to this complex late stage of the evolution, attempts at 
constructing EOB-based waveform models valid up to 
merger [36] and possibly beyond [37, 38] are still in their 
infancy. For NS-BH systems — the focus of this paper 
— one expects relatively large mass ratios, which cause 
complications at both the analytical and numerical level: 
the convergence of the analytical PN approximation is ex¬ 
pected to be slower than for NS-NS systems [39], residual 
eccentricity in the initial data can be appreciable [40-42] 
(but see also [43, 44]), and very different dynamical time 
scales must be tracked by numerical evolution codes [45- 
63]. 

Despite these obstacles, a clear picture of the GW 
emission of NS-BH binaries has emerged over the last few 
years. Most of the GW signal is emitted before the NS is 
tidally disrupted — if this happens at all — and before 
significant thermal effects occur. Furthermore, magnetic 
fields appear to barely affect GW emission [59]. These 
are particularly fortunate circumstances, as they imply 
that an ideal fluid-dynamics treatment with a cold EOS 
and an ideal-gas T-law for the thermal part are appropri¬ 
ate to simulate the dynamical regime that is of interest 
for the GW signal [51]. At least two papers attempted 
a phenomenological description of the GWs emitted by 
NS-BH binaries. 

In the first paper. Lackey et al. [64] developed an an¬ 
alytic representation of the NS-BH IMR waveform cal¬ 
ibrated to 134 numerical waveforms produced by the 
SACRA code [65] with the main goal of assessing the mea¬ 
surability of the NS tidal deformability. 

A subset of these simulations for systems with non¬ 
spinning BHs was then used in work by Pannarale et al. 
[66] (henceforth Paper I) to obtain a phenomenological 
NS-BH IMR waveform amplitude model in the frequency 
domain. This model was, at heart, a “distortion” of the 
PhenomG BH-BH model. Paper I paid particular atten¬ 
tion to the accuracy of the GW spectrum at high fre¬ 


quencies — where the EOS-related phenomenology takes 
place — and to the determination of a cutoff frequency 
in the GW emission. This cutoff frequency is especially 
important in the construction of NS-BH template banks. 
If a BH-BH-like template built to detect a disruptive NS- 
BH coalescence were to be truncated at a frequency that 
is too low with respect to the physical cutoff frequency of 
the source, a loss in recovered signal-to-noise ratio would 
occur. If on the other hand the truncation frequency 
were to be increased in order to counteract this problem, 
it could become too high with respect to the physical 
cutoff frequency of the NS-BH source, and this would 
possibly result in penalizing the template by degrading 
its performance in chi-square tests, which would also be 
detrimental to the detection. 

The goal of the present paper is to extend the work of 
Paper I to NS-BH systems with a non-precessing, spin¬ 
ning BH, using the full set of 134 hybrid waveforms con¬ 
sidered in [64]. The phenomenological model based on 
this catalog allows us to produce the most accurate de¬ 
termination of cutoff frequencies for NS-BH GW signals, 
with relative errors on the cutoff frequency below 10%. 
These errors are well below the errors one would obtain 
using either BH-BH models or the NS-BH model of [64], 
with immediate applications in setting up template banks 
to target these systems. As for Paper I, we adopt a con¬ 
servative approach and focus on the analytical modeling 
of the GW amplitude in the frequency domain, because 
residual eccentricity in our initial data and the short du¬ 
ration of our simulations do not guarantee an accurate 
phasing in the whole parameter space: see Hannam et al. 
[67] for how mass ratio affects the minimum number of 
numerical waveform cycles necessary to ensure an accu¬ 
rate phase and amplitude modelling, and [64] for issues 
in building hybrid waveforms for NS-BH binaries. 

The plan of the paper is as follows. In Secs. H and 
HI we review the basics of the numerical simulations 
and gravitational waveform hybrids, respectively, used to 
build and test the phenomenological model discussed in 
this paper. In Sec. IV we describe the waveform model 
for NS-BH binaries with a spinning BH. In Sec. V we 
compare our model against numerical data. Sec. VI dis¬ 
cusses some important applications of our model, in par¬ 
ticular predictions for the tidal disruption frequency and 
their implications for GW detection and the modeling of 
short gamma-ray bursts (SGRBs). Finally, in Sec. VH 
we summarize our conclusions and point out directions 
for future work. Throughout the paper, unless otherwise 
noted, we use geometrical units (G = c = I). 


II. THE NUMERICAL SIMULATIONS 

Our phenomenological models are calibrated to and 
tested against the gravitational waveforms used in [64]. 
The waveforms are derived by numerical-relativity sim¬ 
ulations performed by the SACRA adaptive-mesh refine¬ 
ment code [65]. The details of the code are described in 
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[47]. Here we only briefly discuss the key differences with 
respect to the simulations performed to derive the wave¬ 
forms used in Paper I. Binaries in quasiequilibrium states 
are prepared as initial conditions for the simulations us¬ 
ing the multidomain spectral method library LORENE [68] . 
In this work we allow the BHs to have nonzero spins 
aligned with the orbital angular momentum of the binary. 
The formulation and numerical methods for computing 
quasiequilibrium configurations are the same as in [69], 
except for the implementation of BH spins [47]. Grav¬ 
itational waveforms are computed from the Weyl scalar 
'I '4 by integrating twice in time using a so-called fixed- 
frequency integration method [70] to filter out unphysical 
low-frequency components (see also [64, 71]). 

We adopt piecewise polytropic EOSs, which mimic 
nuclear-theory-based EOSs with a small number of pa¬ 
rameters [72], to model the NS matter at zero tempera¬ 
ture. Each piecewise polytrope is characterized by poly¬ 
tropic constants Ki and adiabatic indices as 

P{p) = Kip^' for < p < Pi (t = 1 ... , n), (1) 

where p and P are the rest-mass density and pressure, 
respectively. At the critical densities pi we further require 
the pressure to be continuous, i.e. 

^iPi ' ~ ^i+lPi ^ ) ( 2 ) 

and the EOS is thus completely specified by ki, E^, and 
Pi {i = 1 ... , n). In this work, we adopt the same piece- 
wise polytropes that were adopted in Paper I and in [64] . 
More specihcally, n is set to be 2, the parameters {ki, Ei} 
for the low-density crust regions are fixed, and the two 
parameters r 2 and Pfldu = P{p = g/cm^) are sys¬ 

tematically varied to span a plausible range of nuclear- 
matter properties. In the dynamical simulations, thermal 
corrections are added in an ideal-gas-like form in order 
to capture the effect of shock heating [46, 47]. 

The 134 NS-BH simulations used in this paper 
are listed in Table II of [64]. The mass ratio 
Q = Mbh/Mns spans the values {2,3,4,5} and the 
BH dimensionless spin parameter x takes the values 
{-0.5,0,0.25,0.5,0.75}. We adopt x = -0.5 only for 
Q = 2, because the combination of negative x and large 
Q (say > 3) yields small tidal effects during the coa¬ 
lescence. The NS mass Mns is set to 1.35Mq for all 
the runs, with the exception of some {Q,x) = (2,0.75) 
and (2, 0) cases, in which Mns can also take the values 
{I. 2 OM 0 , 1 . 45 M 0 } and I. 45 M 0 , respectively. The EOSs 
are the same 21 models used in Paper I and [64] (see 
Fig. 1 therein for a representation in the piecewise poly¬ 
tropic EOS parameter space) and Fig. 1 in [46] for the 
NS equilibrium sequences yielded by these EOSs). For 
all combinations of Q and x with Mns = 1.35M0, the 
runs are performed adopting EOSs with r 2 = 3.0; addi¬ 
tionally, EOSs with r 2 = 2.4, 2.7, and 3.3 are employed 
for models with (Q, x) = (2,0), (3,0.5), and (5,0.75). 

To build and test our phenomenological NS-BH wave¬ 
form model we use the hybrid waveforms of [64], which 


are also based on the numerical-relativity simulations just 
described (see next section). As in Paper I, we divide 
the datasets into two groups: one to build the waveform 
model and one to test it. In the simulations used to 
build the model the NS mass is Mns = I.35M0, and we 
use the r 2 = 3.0 EOSs denoted by 2H, H, HB, and B 
with log(Pfidu/(dyne/cm^)) = 34.9,34.5,34.4, and 34.3, 
respectively (see, for example. Paper I for this nomencla¬ 
ture): these are 59 datasets. The remaining 75 cases are 
not used to tune the waveform model, but just to test it. 


III. THE HYBRID WAVEFORMS 

In order to build our phenomenological NS-BH 
frequency-domain waveform amplitude model, we must 
first construct accurate IMR waveforms. This is done by 
matching each of the numerical NS-BH waveforms de¬ 
scribed in the previous section — which all begin ~ 10 
GW cycles before merger — to an inspiral waveform 
model, and by then splicing them together. We use the 
PhenomG BH-BH model of [11] as our inspiral waveform 
— a sound approximation, as tidal effects on the am¬ 
plitude are negligible in this stage — and, unless other¬ 
wise noted and in accordance with the conventions of 
[11], all frequencies in this section and in the rest of 
the paper are to be intended as multiplied by the sum 
mo = Mns + AIbh of the two initial masses, i.e. we use 
units in which mo = 1. Similarly, times are to be in¬ 
tended as divided by mo. 

When matching waveforms, a time constant r and 
phase constant (f) are the two free parameters that need 
to be fixed. For a generic waveform h(t), the time and 
phase can be adjusted to produce a shifted waveform 
(t; r, (()) = h{t — . The Fourier-transformed 

waveform^, which can be written in terms of ampli¬ 
tude and phase as h{f) = A(/)e®*^'^\ has a correspond¬ 
ing shifted waveform /i®’"^‘(/; r, 0) = A(/)e*^ 
where r, (()) = 4)(/) -|- 27r/r -|- (j). Because the 

time and phase constants have no impact on the GW 
amplitude in the frequency domain, we do not need to 
calculate them in this work, where we only model the 
amplitude A{f). 

Although the time and phase constants do not impact 
the amplitude, there is still some freedom in construct¬ 
ing the hybrid. We will use the method of Lackey et 
al. [64]. In the time domain, the numerical waveform 
begins with a finite amplitude, leading to the oscillatory 
Gibbs phenomenon that results from Fourier transform¬ 
ing a waveform segment with nonzero starting amplitude. 
We therefore begin by windowing the numerical NS-BH 
waveform with a Hann window over the interval Wi to Wf 


^ We omit the tilde over Fourier-transformed quantities in order 
to keep the notation lighter, as h{t) will no longer appear in the 
rest of the paper. 
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(of width Wf — Wi) defined as 


Won{t) 


1 

2 


1 — cos 


/ Tr[t - 
\wf-wj 


(3) 


We set the start of the window to be the start of the 
numerical waveform at = 0 and use a width of 300 by 
setting Wf = 300, as in Ref. [64], 

We are also free to choose the frequency interval for 
splicing the numerical waveform onto the analytic, inspi¬ 
ral waveform. This interval should be at high enough 
frequencies to exclude the effects of windowing at the be¬ 
ginning of the numerical waveform. It should also exclude 
the small initial eccentricity (eo 0.03) that dies down 
after a few orbits and results from providing the numer¬ 
ical simulation with quasicircular (zero radial velocity) 
initial conditions that ignore the small radial velocity due 
to GW radiation reaction. However, the splicing inter¬ 
val should also be at a frequency that is low enough to 
capture the matter effects, present in the numerical simu¬ 
lations, that take place just before merger. We smoothly 
turn on the numerical waveform and smoothly turn off 
the analytic, inspiral waveform within a splicing window 
Si < f < sf using Hann windows: 


Wosif) 


Wonif) 


1 

2 

1 

2 


1 -f cos 
1 — cos 


7r[/ - 5i] \ 
sf-s, ) 


^[/ - gi] Y 
sj-sv ). 


The amplitude of the hybrid waveform is then 


(4) 

(5) 


^hybrid (/) — 


^bbh(/), 

WoeiD^BEKif) + Won(/)^NR(/), 
^nr(/), 


/ < Si, 

Si< f < Sf, 

f > Sf- 

( 6 ) 


As in Ref. [64], we use a starting frequency of Si = 0.018 
and an ending frequency of Sf = 0.019. 


IV. MODELING SPINNING NEUTRON 
STAR-BLACK HOLE WAVEEORMS 


centered in /o with width d. 

Before discussing the waveform model itself, we must 
introduce two reference GW frequencies: these are /tide 
and /rd ■ The former is the GW frequency at the onset of 
the NS tidal disruption, while the latter is the dominant 
(£ = m = 2, n = 0) ringdown frequency of the remnant 
BH. The BH remnant dominant ringdown frequency /rd 
depends on the the mass Mf and spin parameter xi of the 
BH remnant of the NS-BH merger. These are calculated 
according to the model discussed in [73, 74], while the fit¬ 
ting formulas that relate xt and Mf to /rd are provided 
in [75]. To calculate the parameters of the BH remnant 
given the initial parameters of the binary, we follow the 
model detailed in [73, 74]. To compute /tide, on the other 
hand, one must first determine a coefficient ^tide that pro¬ 
vides a relativistic correction to the standard Newtonian 
estimate of the orbital radius at mass-shedding [76]. This 
coefficient can be found by solving the equation 

-^NsCtide _ 3[^tide ~ ^M^tide + /g\ 

Ctide - 3/X^tide + 2X\/M^6ide ’ 

where fi = Mbh/Rns = QC, with C = Mns/Rns de¬ 
noting the NS compactness. The orbital radius at mass- 
shedding is now given by 


Uide = ,CtideRNs(l — 2C) . (9) 


The tidal frequency /tide then reads 

/tide = ±———- ,.o t ^ (Id) 

7r(xfMf -t 

where upper/lower signs hold for prograde/retrograde or¬ 
bits. So far — aside from the inclusion of y-dependent 
terms in Eqs. (8)-(10) and in obtaining Mf and Xf from 
the model of [73, 74] — nothing differs from the approach 
laid out in Paper I. We would like to note that in the pro¬ 
cess of calculating Xf and Mf one must also determine 
another quantity that plays a role in the gravitational 
waveform model. This is the mass of the torus that may 
remain around the BH at late times, Mb,torus, modelled 
using the fitting formula [76] 


In this section we provide a detailed description of 
our phenomenological model for the frequency-domain 
GW amplitude of nonprecessing NS-BH binaries with a 
spinning BH component. This new model generalizes 
the model presented in Paper I for nonspinning binaries; 
throughout the discussion we will point out differences 
with respect to the formulation reported in Paper I. 

In accordance with the simulations and hybrid wave¬ 
forms at our disposal, we set the BH spin vector to be 
aligned to the orbital angular momentum of the binary. 
Additionally, we use the notation ^(/) for the win¬ 
dowing functions 


w 


± 

fo,d 


if) 


1 

2 


I ± tanh 




(7) 


,torus _ 0.296ftide — 0-17Irisco 
Mb,NS Rns 

where Mb ,ns is the rest-mass of the NS in isolation and 
risco is the radius of the innermost stable circular orbit 
(ISGO) of the initial BH in isolation [77]: 


Usco — [3 -I- .^2 T '\/(3 — ^i)(3 Zi 2 Z 2 )], 

Zf = I + (1 - x 2)1/3 + ^)l/3 + (1 _ ^)l/3 


Zn = 


( 12 ) 


As in Paper I, we write the amplitude Apben(/) of the 
frequency-domain signal hphen(/) = Aphen(/)e*^''*'“(7) 


























5 


as a sum of three terms: 

+ (I!*) 

where Apn(/) is the inspiral contribution, based on the 
stationary-phase approximation and obtained by combin¬ 
ing a 3PN-order time-domain expansion of the amplitude 
and the TaylorT4 description for the phase (see [11] for 
further details); the second term models the premerger 
and merger (strong-held) modihcations to the PN inspi¬ 
ral (the 7i coefficient is provided, once more, in [11]); and 
-^RD is the ringdown amplitude. This is modelled via a 
Lorentzian £(/, /o, a) = cr‘^/[{f - fof + cr^/4]: 

AMf) = etide^i/:(/,/RD(Xf,Mf),5'/RD/Q(Xf))/-"/®, 

(14) 

where^ (5i is the ringdown amplitude htted to BH-BH hy¬ 
brid waveform data in [11], CtMe is a NS-BH correction 
discussed later on in this section, and is a fudge factor 
which accounts for errors in the model used to compute 
Xf, as this spin parameter is in turn used to determine 
the quality factor Q of the BH remnant (once again via 
the htting formulas of [75]). We note that the PhenomC 
model also uses a fudge factor, S 2 , but Xf is determined 
using the formulas in [78]. Therefore, rather than cor¬ 
recting the PhenomC parameter 62 as we did in Paper 
I, in this paper we introduce a S 2 parameter in order to 
disentangle PhenomC and our model more clearly. The 
remaining elements of Eq. (13) to be discussed are the 
windowing functions. As for the PhenomC model, d is 
set to 0.015, but, just as in Paper I, we allow for a cor¬ 
rection tTtide to the width, we do not necessarily tie the 
central windowing frequencies to the BH remnant ring- 
down frequency, and we do not always fix the central 
windowing frequencies of the first two terms in Aphen(/) 
— Cins/o and /o, respectively — to the same value. More 
details on how crtidei eins and /o are determined are given 
further on in this section. 

To summarize, in order to build Aphen(/) given Mbh, 
X, Mns and an EOS (which determines the NS radius 
i?Ns and its baryonic rest-mass Mb,Ns); one must be¬ 
gin by computing: (1) 71 and (5i according to [11], (2) 
Afb,torus using Eqs. (11)-(12), (3) Mf and Xf following the 
model reported in [73, 74], (4) fwiMf, Xf) and Q(xf) via 
the fits of [75], and (5) /tide, following Eqs. (8)-(9). At 
this point the model splits into four cases, depending on 
the values of /tide, /rd, and Mb,torus- These cases reflect 
the different phenomenology of NS-BH binary mergers 
observed in the 59 simulations used to build the model 
(see Sec. H): “disruptive,” “nondisruptive,” and “mildly 
disruptive” with and without a torus remnant. 


^ The third argument here corrects a typo in Eq. (3) of Paper I. 


Before laying out the necessary details about the four 
alternative waveforms in the remaining subsections, we 
will briefly explain how the four phenomenological mod¬ 
els were obtained. Each of the hybrid waveform ampli¬ 
tudes corresponding to nondisruptive mergers — out of 
the 59 hybrids used to build our model — was fitted with 
the ansatz in Eq. (15) below, leaving etide, Ctide, and $2 as 
free coefficients. The values of the free coefficients were 
then themselves fitted as detailed in Eqs. (16)-(20). The 
same procedure was followed for the disruptive mergers, 
where, this time, the waveform amplitude ansatz is given 
in Eq. (21) and the values of the free coefficients Cms and 
CTtide are fitted with Eqs. (22)-(25). All equations involved 
in this process were inspired by the nonspinning study re¬ 
ported in Paper I. For the two mildly disruptive cases, 
we adopted a strategy that is similar to the one discussed 
in Paper I: we picked and combined ingredients from dis¬ 
ruptive and nondisruptive phenomenological waveforms, 
without having to perform additional fits. 


1. Nondisruptive Mergers 


If /tide > /rd and Mb,torus = 0, the merger is “nondis¬ 
ruptive.” Notice that the first condition differs slightly 
from the one used in the nonspinning waveform model re¬ 
ported in Paper I. Additionally, we had not explicitly re¬ 
quired the second condition in Paper I, because we never 
encountered nonspinning cases for which the analytical 
models would predict Mb, torus > 0 and /tide > /rd- For 
nondisruptive mergers the binary can complete its full 
inspiral stage, and therefore we set eins = 1 in Eq. (13). 
Moreover, the BH remnant ringdown is excited, so the 
windowing functions in Eq. (13) can be centered around 
/o = /rd^- Aphen(/) thus takes the form 


Aphen(/) = Apn(/)w 


/RD5<i+<^tide 

+ 1-257i/ / 

+ '^RD(/)<...+.«ae’ 


(15) 


where Ard is given by Eq. (14). 

The parameters entering Eqs. (14) and (15) are given 
by: 


etide = 2w+ di (2 ^nd) , (16) 


^ Nondisruptive NS-BH mergers essentially behave as BH-BH 
mergers. In PhenomC waveforms, /o is set to O.QS/rdj this 
quantity in turn depends on the mass of the BH remnant, which 
is set to the sum of the individual masses of the binary compo¬ 
nents, rriQ. In our model the mass of the BH remnant is instead 
determined with the approach described in [73]. More specifi¬ 
cally, Mf is always smaller than mo, the final mass choice made 
in PhenomC. This increases /rd? hence the need for the extra 
0.99 in our definition of /q. 
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Q=5., Xbh=0.5, /Wns=1.35Mq, EOS H 



Q=2.,Xbh=0.75, /Wns=1.35Mq, EOS HB 



Q=2., Xbh=-0.5, Mns=1.35/Wq, EOS B 



FIG. 1. Examples of nondisruptive, disruptive, and mildly 
disruptive merger waveforms (from top to bottom). The hy¬ 
brid waveform (in gray) is compared to the PhenomC model 
(dot-dashed blue), to the nonspinning NS-BH model of Paper 
I (dotted red), to the NS-BH model of [64] (dashed green), 
and to the model developed in this paper (long-dashed or¬ 
ange). The short vertical lines denote /rd and /tide in blue 
and dotted red, respectively. 


where xi = —0.0796251, di = 0.0801192, and 


a^ND 


/tide /rd 

/rd 


0.00508451x; 


0.571505C 


(17) 


O’tide — ^'^x2,d2(^ND) ) 

where X 2 = —0.206465, ^2 = 0.226844 and 


a^ND 


f -0.657424C 

\ /rd / 

0.0259977X; 


and finally 


^2 = 


/tide /rd 

/rd 


(18) 


(19) 


( 20 ) 


where A = 1.62496, X 3 = 0.0188092 and ds = 0.338737. 
Eq. (16) slowly suppresses the ringdown of the BH rem¬ 
nant as the merger becomes less and less similar to the 
BH-BH case. We notice that, at variance with the non¬ 
spinning formulation of Paper I, we had to introduce two 
different independent variables for the two waveform pa¬ 
rameters Ctide and (Ttide; i-G. Xnd 7 ^ ai^Q. Furthermore, 
the two independent variables now contain terms linear 
in X- An example of nondisruptive merger spectrum is 
shown in the first panel in Figure 1. Here, the hybrid 
waveform (in gray) is compared to the PhenomC model 
(dot-dashed blue), to the nonspinning NS-BH model of 
Paper I (dotted red), to the NS-BH model of [64] (dashed 
green), and to the model developed in this paper (long- 
dashed orange). 


Disruptive Mergers 

If /tide < /rd and Mbyorus > 0 , then the merger is 
“disruptive,” the NS material is scattered around the 
BH, and the ringdown contribution to Eq. (13) vanishes, 
i.e. Ctide = 0 in Eq. (14). As in the case of nondisruptive 
mergers, the definition of this class is slightly different 
from its corresponding definition in Paper I. For disrup¬ 
tive mergers the waveform model reads: 

Aphe„(/) = ApN(/)^C-_^^/^,,^,d+aRde 

+ l-257i/^/^«^7,,..d+<R.,.> (21) 

which is equivalent to Eq. (13) with /o = /tide and etide = 
0. The remaining parameters to be prescribed are eins 
and (Ttide- These are given by: 

^ins — 


Oi -I- &iXd 


( 22 ) 
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where ai = 1.29971, bi = —1.61724 and 


xd = 


,torus 

A4^b.NS 


+ 0.424912C + 0.363604yi; 


- 0.0605591X, 


(23) 


V = M^sM-bh/'itT'o being the symmetric mass ratio; and 
CTtide = 02 + 62X0 , ( 24 ) 


where 02 = 0.137722, 62 = —0.293237 and 

ccd = - 0.132754C + 0.576669^11 

tWb.NS 

- 0.0603749X - 0.0601185x2 

- 0.0729134x^ . (25) 

As in the case of nondisruptive mergers, two distinct def¬ 
initions of the independent variable are used for the two 
waveform parameters, which was not the case in Paper 
I. Incidentally, this allows us to simplify the functional 
form of the parameter (Jtide with respect to Paper I. An 
example of disruptive merger spectrum is shown in the 
middle panel of Figure 1. Notice that the nonspinning 
model of Paper I is completely inaccurate when high spin 
and strong tidal effects come into play, as the model did 
not account for the combination of the two. 


3. Mildly Disruptive Mergers with no Torus Remnant 

If /tide < /rd and Mb,torus = 0 , then the merger is 
“mildly disruptive,” but no remnant torus is formed dur¬ 
ing the coalescence. The behavior of the GW amplitude 
for this class of mergers is captured by using the disrup¬ 
tive merger waveform model (Sec. IV 2) and modifying 
two of its features. The first modification concerns the 
central frequencies of the windowing functions entering 
Eq. (21). In the present case we set 

Aphenif) = ApN(/)w7^,d+atide 

+ l-257i/^/^a;7^.d+.,,ue> (26) 

where /i = (I - (5 “^)/rd + Q“^ei„s/tide and /2 = 
(1 — Q“^)/rd + Q~^ftide (our definition of the binary 
mass ratio is such that Q > 1). These two frequencies are 
introduced in order to obtain a smooth transition from 
nondisruptive to disruptive merger waveforms, which use 
/rd, and Cins/tide and /tide as windowing central frequen¬ 
cies, respectively. This constitutes a major improvement 
with respect to Paper I, where the low number of avail¬ 
able nonspinning mildly disruptive simulations had not 
allowed us to go into such fine modelling details. Addi¬ 
tionally, the scaling with Q is such that the model can 
accurately reproduce the data for mildly disruptive merg¬ 
ers with no torus remnant. 

The second modification with respect to the disruptive 
merger waveform model in Eq. (21) is in the calculation 


of (Ttide- For mildly disruptive mergers with no torus 
remnant, this calculation is performed by averaging the 
disruptive and the nondisruptive prescriptions: i.e., we 
evaluate Eqs. (18) and (24) and divide their sum by two. 
As for /i and / 2 , this modification allows for a smooth 
transition from disruptive to nondisuptive merger wave¬ 
forms, and it is a notable improvement with respect to 
Paper I. 


4- Mildly Disruptive Mergers with a Torus Remnant 


If /tide > /rd and Mr, torus > 0 , then the merger is 
“mildly disruptive” and a small remnant torus is pro¬ 
duced. Mathematically, this class of mergers arises from 
a shortcoming of the approximations that lead to /tide, 
Mb,torus, and, possibly, /rd: one would not expect a rem¬ 
nant torus to be formed if the system cannot reach the 
onset of tidal disruption. Physically, what the combina¬ 
tion /tide > /rd and Mb,torus > 0 suggests (and what 
the hybrid waveforms that fall in these categories indeed 
show) is that some mildly disruptive mergers can achieve 
both a quasinormal mode (QNM) excitation of the BH 
remnant and the formation of a small remnant accretion 
torus^. From a phenomenological point of view this is 
to be expected, as the outer layers of the NS may be 
stripped to form the torus, while the core of the NS (or a 
fraction of it) may accrete onto the BH coherently enough 
to trigger the BH QNM oscillations (with regards to this 
topic, see Fig. 17 in [47] and our discussion of Fig. 6 be¬ 
low) . We wish to remark that this level of sophistication 
in cataloguing mildly disruptive mergers was not possible 
in Paper I, and is dictated by the fact that in the present 
work we have a larger catalog of numerical-relativity data 
to reproduce. 

For mergers with fade > /rd and Mb,torus > 0, the 
waveform amplitude model reads 


Aphen(/) = Apn(/)ui 


€ins/RDiCi+»^tide 

+ ^RD(/)^L.<^+-.na.’ 


(27) 


where eins is set according to Eq. (22), as for disruptive 
mergers, and (Jtide is given by Eq. (18), as for nondis¬ 
ruptive mergers. Notice that Eq. (27) resembles Eq. (15), 
with the exception of the use of cins in setting the first two 
windowing central frequencies: this allows for a smooth 
transition between nondisruptive and disruptive wave¬ 
form models whenever the NS is tidally disrupted, but 
the excitation of the BH remnant ringdown takes place 
and no torus remnant is expected to survive the merger. 
An example of disruptive merger spectrum is shown in 
the last panel in Figure 1; in this specific case, Eq. (27) 
is used to model the waveform amplitude. 


^ Mb torus was found to be < O.OIMq for these tuning cases: this 
is clearly below the precision of Eq. (11). 







FIG. 2. Relative error distribution over the 134 NS-BH merger simulations for /cut (left panel) and h{fcut) (right panel). All 
134 results from the hybrid waveforms are compared to the predictions obtained from different models: the waveform model 
reported in this paper (“PhMix,” orange), PhenomC (“PhC,” blue), the waveform model of [64] (“Lackey et ah,” green), and 
the simple estimate given by min (/rd,/ tide) (purple). As discussed in the text this last proxy does not provide a waveform 
amplitude model, so it does not appear in the panel on the right. 


V. TESTING THE WAVEFORM MODEL 

Once the phenomenological NS-BH waveform ampli¬ 
tude model is formulated and calibrated via fits to hy¬ 
brid waveforms built upon numerical-simulation data, its 
accuracy must be tested and demonstrated. 

The qualitative behavior of the spinning waveform 
model is the same as for the nonspinning model formu¬ 
lated in Paper I. For the sake of brevity, we do not show 
all 134 GW spectra. We perform, instead, a quantitative 
assessment of the accuracy of the model in terms of two 
characteristic frequencies of each GW spectrum — /viax 
and /cut — and the amplitude of each spectrum at these 
two frequencies. For the cutoff frequency /cut, we adopt 
the same, general definition introduced in Paper I, which 
also involves introducing and determining /mux- This 
definition has the advantage of allowing for a straight¬ 
forward comparison among GW spectra originating from 
different models and/or calculations for the same binary, 
and for consistent comparisons among binaries with dif¬ 
ferent physical parameters. The definition is as follows: 
/Max is the frequency such that f^h{f) has a maximum, 
and /cut is the frequency (greater than /Max) at which 
the amplitude drops by one e-fold, i.e., 

efcuMfcut) = /Max/i(/Max) ■ (28) 

We stress once more that this definition of the cutoff 
frequency is independent of the details of the waveform, 
and it works for any h{f) (given in either analytical or 
numerical form). 

Figures 2 and 3 show the distributions of the relative 
errors (e) on /cut and /.(/cut), and /Max and /(/Max), 
respectively. For each of the 134 hybrid waveforms, we 
determine these four quantities and compare them to the 
values predicted for them by several models. These are: 
(1) the waveform model reported in this paper (labelled 


“PhMix” in the Figures), (2) PhenomC (“PhC”), (3) the 
waveform model of [64] (“Lackey et al.,”), and (4) the 
simple prescription min(/FiD, /tide), which, of course, can¬ 
not be used to predict /(/cut) and /(/Max), but just as 
a tentative proxy for /cut and /Max- While the low fre¬ 
quency, inspiral regime is described by construction in 
the same way by all models, our new phenomenologi¬ 
cal model clearly introduces a considerable improvement 
in terms of accurately predicting all four high-frequency 
features of GWs emitted by spinning NS-BH coalesc¬ 
ing binaries. Both the mean and maximum values of 
the e(/cut) and e(/Max) distributions are considerably re¬ 
duced when using our model. The maximum relative 
error on the cutoff frequency is of order ^ 10%, to be 
compared with ^ 60% for the model of [64] and the 
proxy min(/RD, /tide), and to even higher relative er¬ 
rors for PhenomC, which was not designed for NS-BH 
binaries. The simple prescription min(/RD, /tide) is a 
better proxy for /cut than it is for /mux- Tbe advan¬ 
tages of our model are equally striking when considering 
e(/(/cut)) and e(/(/Max))- Both distributions are peaked 
at < 0.05 and fall off rapidly for the phenomenological 
NS-BH model, whereas this is not the case for PhenomC 
and the model of [64]. In conclusion, our approach is 
more accurate than all of these existing alternatives in 
modelling the GW spectra of NS-BH binaries. 

In Figure 4 we separate the two contributions to the 
e(/cut) distribution shown in Figure 2: in the left panel 
we show the distribution of the relative error on the cutoff 
frequency for the 59 binaries use to calibrate the model, 
while in the right panel we report the same distribution 
for the 75 binaries that were only used to test the model. 
A comparison with the left panel of Figure 2 shows that 
the performance of our model is not dominated by results 
for “calibration binaries.” Analogous behaviors were ob¬ 
tained for e(/Max), e(/(/cut)) and e(/(/Max))- 
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FIG. 3. Relative error distribution over the 134 NS-BH hybrid waveforms for /Max (left panel) and /i(/Max) (right panel). 
Colors and labelling are the same as in Figure 2. 



0.0001 0.001 


0.01 0.1 
e(fcut) 


1 . 10 . 



FIG. 4. Relative error distribution on /cut (as in Figure 2), but now the left panel refers to the 59 “calibration binaries” used 
to build the model, while the right panel refers to the 75 “test binaries” that were not used to build the model. 


VI. APPLICATIONS 
A. The Cutoff Frequency 

The panels in Figure 5 show the cutoff frequency /cut 
of nonspinning NS-BH binaries as a function of the NS 
mass Mns and the binary mass ratio Q for the four dif¬ 
ferent EOSs (2H, H, HB and B) used to compute the 
hybrid waveforms on which we built our phenomenolog¬ 
ical GW amplitude model. This figure is an updated 
version of Figure 9 in Paper I. The contour line values 
are reported in Hz, with a 250 Hz spacing. The two 
dashed lines in each panel separate disruptive mergers 
(bottom-left region), non-disruptive mergers (top-right 
region), and mildly disruptive mergers (middle region in 
between the two lines). We remind the reader that dis¬ 
ruptive mergers correspond to the conditions /tide < /rd 
and Mb,torus > 0; nondisruptive mergers are such that 
/tide > /rd and Mb,torus = 0; and mildly disruptive 
mergers do not fall into either of the previous categories. 
To help the reader, the three regions are explicitly in¬ 


dicated in the top-left panel of the figure. In all pan¬ 
els, a green line marks the boundary between binaries 
with Mb,torus > 0 and those with Mb,torus = 0. Notice 
how these green lines overlap with the lower dashed lines: 
this indicates that the mildly disruptive mergers in these 
panels are such that the NS is tidally disrupted, but no 
remnant torus is formed. The main difference between 
these plots and those in Paper I is that the behavior of 
/cut across the dashed transition lines is now smoother, 
and that the contours are continuous. This improvement 
over the original model is due to the larger data set used 
here, which allowed us to better tune our model in the 
mildly disruptive region. 

The spin dependence of /cut is illustrated in Figure 6, 
where the different panels refer to initial BH spin parame¬ 
ters X S {—0.5, —0.25,0,0.25,0.5, 0.75} and we show two 
extreme cases for the EOS (2H and B). By comparing the 
EOS 2H and the EOS B panels, we therefore get an idea 
of the span of possible cutoff frequencies /cut at a given 
BH spin value. As expected, the main differences oc¬ 
cur in the disruptive merger regions, as this is where the 
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EOS 2H, x=0 EOS H, x^O 




FIG. 5. The cutoff frequency /cut, as defined in Eq. (28), computed with our NS-BH GW amplitude model. We set the BH 
spin parameter x to zero and consider the EOSs 2H, H, HB, B. The contour lines report /cut in Hz and have a spacing of 
250 Hz. The thick, green, continuous line is the location where Mb,torus vanishes, that is, where the left hand side of Eq. (11) 
is zero. The dashed lines in each panel divide the plane in three regions, explicitly labelled only in the top-right panel to avoid 
overcrowding the plots. These are: a top-right region in which the NS-BH coalescences are nondisruptive, i.e. /tide > /rd and 
Mb,torus = 0; a bottom-left one in which they are disruptive, i.e. /tide < /rd and Mb,torus > 0; and a middle region in which 
mildly disruptive coalescences occur, i.e. /tide < /rd and Mb,torus = 0, or /tide > /rd and Mb,torus > 0. This Figure should be 
compared to Fig.9 of Paper I: notice how the contour lines and the transitions between different regimes are smoother with the 
new model we formulate in this paper. 


NS EOS impacts the dynamical evolution of the binary. The relative size of the disruptive region grows with y. 
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EOS 2H, x=-0.5 


EOS 2H, x=-0.25 


EOS 2H, x=0 











FIG. 6. Same as Fig. 5, but focusing only on the “extreme” EOSs 2H and B. The BH spins range between —0.5 and 0.75 in 
steps of 0.25, as indicated in each panel. 
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because larger spins increase the likelihood of tidal dis¬ 
ruption. It is interesting to notice that above x — 0-5 
the green torus mass boundary lines no longer track the 
lower dashed lines, but the upper ones (cf. the panels 
with X = 0.75, and the EOS B panel with y = 0.5). 
This means that a small remnant torus is likely to be 
formed for mildly disruptive mergers in these cases. It 
also demonstrates that there is indeed a need to split the 
phenomenological GW model into four subcases, as we do 
in this paper, at least until better analytical predictions 
are available for Mb^torus) /tide, and /rd, if one wants to 
keep using these as tools to build the waveform model. 


B. Phenomenological Fits 

As shown in a companion paper, the results presented 
above can be used to construct simple phenomenological 
formulas to determine (1) whether an NS-BH binary is 
disruptive, mildly disruptive, or nondisruptive, and (2) 
the GW cutoff frequency for disruptive mergers, due to 
the NS tidal disruption. Details are given in [79], and 
here we only review the main results. 


1. Disruption Criterion 

The contours that separate NS-BH binaries with a 
disruptive fate from those with a mildly disruptive or 
nondisruptive fate in Figs. 5 and 6 may be fitted in sev¬ 
eral ways as a function of the physical parameters of the 
binary. The critical binary mass ratio Qd = Qd(C,x) 
below which mergers are disruptive is an approximately 
“universal” (i.e. EOS-independent) function that is well 
fitted by: 

3 

Qd= ’ (29) 

i,j=0 

i+j<3 

where the a^’s are constants (see [79] for details). 

In building template banks for GW detection and for 
other applications, it may be useful to know Qd as a 
function of the NS mass for NSs with large radii, as this 
is the case in which the GW emission from NS-BH sys¬ 
tems differs the most from BH-BH binaries. A fit of the 
disruptive boundary for the 2H EOS yields 

3 

Qd= ^ (30) 

ij=0 

i+j<3 

where Mns = Af^s/Af© and the coefficients take the val¬ 
ues reported in Table I in units of G = c = toq = 1- We 
remark again, that this fit gives a lower limit on Qd, as 
it was obtained using an exceptionally stiff EOS. 


2. Cuto ff Frequency Fitting Formula 

When tidal disruption occurs, as determined via 
Eq. (29) , our phenomenological model allows us to deter¬ 
mine a formula that provides the GW cutoff frequency 
analytically, as follows. We consider the 2H EOS, gen¬ 
erate a set of 10^ random disruptive mergers, compute 
/cut for each NS-BH binary according to the definition 
in Eq. (28), and finally fit the data thus obtained. In or¬ 
der to select disruptive mergers, we randomly sample the 
parameter space in the ranges Mns/AT© G [1.2,2.83] and 
Q G [2,10], X G [—0.5,0.75]; we verify whether the sam¬ 
pled point corresponds to a disruptive binary, as defined 
just above Eq. (21), and keep the point if it does. The 
whole process is repeated until we have the desired set 
of 10^ disruptive binaries. While the maximum NS mass 
for the 2H EOS is ~ 2.83M©, the maximum NS mass in 
our sample of disruptive NS-BH mergers is ~ 2.28M©. 
The resulting mass interval A/ns/M© G [1.2,2.28] corre¬ 
sponds to the compactness interval 0.117 < C < 0.221. 
With this set of disruptive cutoff frequency data in hand, 
we fit /cut in terms of the NS-BH binary parameters using 
the ansatz 

3 

/cut - E , (31) 

ij=0 

i+i+k<3 

where the fitting coefficients can be found in [79]. The 
relative error distribution for this ht with respect to the 
original data points is shown in the left panel of Figure 
7. Notice that the peak of the distribution is below the 
percent level: the relative error for 68%, 95%, and 99.7% 
of the points is 0.47%, 1.5%, and 4.9%, respectively. 

As a consistency check for this fitting formula, we draw 
a separate sample of 10'^ disruptive mergers, compute the 
GW amplitude cutoff frequency for each binary, and de¬ 
termine the relative errors of the fit just discussed against 
these “test” binaries. This time we use EOS B to con¬ 
struct our “test” sample and we lower the maximum al¬ 
lowed NS mass to 2M©, as this is approximately the max¬ 
imum Mns for this EOS. The compactness now ranges 
between ~ 0.161 and ~ 0.225. The result of this test is 
reported in the right panel of Figure 7, where we show the 
relative error distribution for the fit in Eq. (31) with re¬ 
spect to the “test” set of binary mergers populated using 
EOS B. Remarkably, the maximum relative error is 2.2%, 
97.6% of the points have a relative error that is smaller 
than 1%, and the peak of the distribution is once again 
below the percent level. The relative errors for the EOS 
B “test” set are even better than for the “calibration” set 
of EOS 2H, because EOS B covers a narrower range in 
compactness relative to EOS 2H. Furthermore, the fit of 
Eq. (31) is effectively EOS-independent (or “universal”), 
at least within the parameter space region in which our 
model was calibrated. 

For GW data analysis purposes, we also performed a 
fit of the 2H EOS /cut data in terms of the NS mass. 
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TABLE I. Values of the coefficients of the fits discussed in Eqs. (30) and (32). The number below each coefficient symbol must 
be multiplied by the power of ten in square brackets on the right-hand side of the coefficient symbol. The Qijk’s are reported 
in G = c = mo = 1 units. 


boo [10^] 

feio [10^] 

boi [10^] 

520 [10^] 

fell [lO’j 

bo2 [lO’] 

bao [10“] 

521 [10“] 

5i2 [10“] 

5o3 [10 ’] 

4.14730 

-5.70783 

2.57882 

2.91134 

-2.44263 

1.04225 

-5.26102 

6.28215 

-5.13944 

3.99706 

gooo [10“^] 

9100 [10“^] 

9010 [10“^] 

ffooi [10“'^] 

3200 [10“'^] 

3020 [10“®] 

3002 [10“^] 

3110 [10“"] 

3101 [10“^] 

3011 [10“'‘] 

9.17677 

-1.39031 

-2.61399 

2.43286 

5.46375 

1.74490 

4.16418 

2.31878 

-7.49673 

5.65265 

9300 [lO”"^] 

5Q30 [10“'’] 

9003 [10“’’] 

3210 [10“'’] 

3120 [10“'‘] 

3201 [10“"'] 

3102 [10“’'] 

3021 [10“'‘] 

3012 [10“'’] 

3111 [10“’] 

-1.28721 

-6.88240 

-2.32757 

3.64301 

-4.05234 

-8.98986 

-1.87475 

5.50808 

-5.78858 

-8.55090 




FIG. 7. Distribution of relative error on /cut obtained by comparing the values given by the fit in Eq. (31) to the lO"^ EOS 2H 
datapoints, used to produce the fit (left) and to 10"^ EOS B datapoints not involved in producing the fit (right). 


rather than its compactness: 

3 

/cut= 51 (32) 

i,j=0 

i+j+k<3 

where the coefficients gijk are listed in Table I. The re¬ 
sulting error distribution is similar to the one shown in 
the left panel of Figure 7, but now the relative errors 
with respect to the EOS B “test” data are much higher 
than those in the right panel of Fig. 7 and range from 
54% to 65%. In other words, the mass fit in Eq. (32) is 
not EOS-universal. The fit is still useful, as it provides 
a lower limit to the GW amplitude cutoff frequency/cut 
as a function of Mns, Q and y. 


VII. CONCLUSIONS 

In this paper we have extended the phenomenological 
gravitational waveform amplitude model of Paper I to 
NS-BH binaries in which the BH has spin either aligned 
or antialigned with the orbital angular momentum. We 
have classified binaries in four broad categories: “disrup¬ 
tive,” “nondisruptive,” and “mildly disruptive” with and 
without a torus remnant. 


The phenomenological tools developed in this pa¬ 
per can be used in various contexts, as detailed else¬ 
where [79]. Our predictions for the GW amplitude and 
for the cutoff frequency may be used to improve the tem¬ 
plate banks currently used in NS-BH searches, and they 
could also be exploited to build new NS-BH phenomeno¬ 
logical GW phase models and EOB models. More accu¬ 
rate gravitational waveforms improve our chances of de¬ 
tecting GW signals and of extracting information from 
them. Eor example they could provide better constraints 
on the NS EOS [80] and possibly even on the underlying 
theory of gravity [81]. Furthermore, our work allows us 
to pin-point binaries in which tidal effects are relevant. 
These cause the GW signal to deviate significantly from 
a BH-BH one and, possibly, lead to the emission of elec¬ 
tromagnetic radiation. The latter may either be in the 
form of a relativistic jet launched by a hot, massive (> 
O.OIMq) disk produced by the tidal disruption of the NS 
(a scenario that could explain the duration, energetics, 
and estimated event rates of SGRBs [82-84]), or in the 
form of isotropical radiation, i.e. macronovae/kilonovse, 
powered by decay heat of unstable r-process elements 
and by nonthermal radiation from electrons accelerated 
at blast waves between the merger ejecta and the inter¬ 
stellar medium [49, 55, 85-89]. For these reasons, our 
model can have an impact on multimessenger searches 
targeting GW, electromagnetic, and neutrino radiation. 
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as well as important applications in the interpretation of 
future multimessenger observations. 

All our predictions are clearly affected by systemat- 
ics in the initial data for the numerical simulations we 
used, in the numerical evolutions, in the phenomenolog¬ 
ical model itself and the tools it relies on [73, 76], and in 
the htting procedures. We expect these errors to increase 
when the model is extrapolated beyond the parameter 
space in which it was tuned. Future work should extend 
and improve our model in order to include not only the 
GW frequency domain amplitude, but also its phase. It 
should also consider larger values of the BH spin, nonzero 
NS spins and (most importantly) precession effects that 
occur when the BH spin is not aligned with the orbital 
angular momentum (see e.g. [48, 62]). Further, any im¬ 
provements in the underlying phenomenological BH-BH 
model can and should be included in our framework for 
NS-BH systems. In particular, the recent “PhenomD” 
BH-BH model [12, 13] is calibrated to hybrid effective- 
one-body waveforms that use numerical-relativity simu¬ 
lations with mass ratios up to 1:18, and BH dimension¬ 


less spin parameters up to ~ 0.85 (0.98 in the equal-mass 
case). This model is an improvement with respect to the 
“PhenomC” model, which was calibrated up to mass ra¬ 
tios of 1:4, and resolves the technical limitations that may 
be encountered when using it for BH-BH systems with 
mass ratio higher than 1:20 and |x| > 0.9. We plan to 
address all of these issues in the near future. 
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